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Abstract 

Second-order equations in terms of auxiliary variables similar to 
potential and stream functions are obtained by applying a weighted least 
squares formulation to a first-order system. The additional boundary 
conditions which are necessary to solve the higher order equations are 
determined and numerical results are presented for the Cauchy— Riemann 
equations. 
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INTRODUCTION 


Least squares methods are often used to solve systems of first-order 
equations. Classical discretization and iterative procedures can then be 
applied to the resulting system of second-order equations. Simply 
differentiating the equations, however, may lead to difficulties particularly 
if the nonhomogeneous terms are not regular. In this paper the least squares 
method is modified and a system of second-order equations is obtained without 
the need to differentiate the nonhomogeneous terms. Numerical examples are 
presented and the advantages of the present formulation are discussed. 

In Section 2 the classical least squares method is reviewed. The modified 
least squares method is introduced in Section 3. Section 4 gives some 
numerical details on the application of the two methods and in Section 5 we 
summarize our main results. 


2. LEAST SQUARES METHOD FOR A SYSTEM OF FIRST-ORDER EQUATIONS 

For illustrative purposes we consider the equations which describe the 
flow in a straight channel with a circular arc airfoil mounted on its lower 
wall. The governing equations are those of continuity and vorticity: 

(pu) + (pv) = s, (1) 

a y 

u - v = 

y x 


- 0 ), 


( 2 ) 
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where p is the density, and u and v are the components of velocity in 
the x and y directions, respectively. The nonhomogeneous terms s and 
to represent some given source and vorticity distributions in the field. The 
linearized boundary conditions associated with this problem are 

u = Uq along x = 0 and x = A, (3) 

v = Vq along y = 0 and y = h, (4) 

where A and h are respectively the length and height of the channel. We 
shall refer to the boundary value problem (1), (2) with boundary conditions 

(3), (4) as Problem I. This problem is depicted in Figure 1. 



u 


0 


Figure 1. Pictorial Representation of Problem I 




-3- 


Let ft denote the domain {(x,y) :0<x<£,0<y<h} and 3ft its 
boundary. Applying Gauss' theorem to the divergence form of the first 
equation yields 

// V* (p^)dA = }> p(q»ti)ds, (5) 

ft 3ft 

where _q = (u,v) and _n is the unit outward drawn normal to the boundary. 

Alternatively, if the tangential velocity is specified along the boundary 
a different compatability relation should be satisfied, in terms of the 
vorticity. Suppose that, in this case, the boundary conditions are given by 


u = g, (x) along y = 0 

( 6 ) 

u = g 2 ^ x ) along y = h 

v = Vq = 0 along x = 0 and x = £ . (7) 

We shall refer to the boundary value problem (1), (2) with boundary conditions 
(6), (7) as Problem II. This problem is depicted in Figure 2. 
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u = g x (x) 


Figure 2. Pictorial Representation of Problem II 


According to Stokes' theorem, the vorticity and the circulation are 
related as follows 

//a)dA=/)q^*^sds, (8) 

ft 9ft 

where £ is a unit vector tangent to the boundary 9ft and the line integral 
is taken in the anti-clockwise direction. 

An application of the least squares method to Problem I leads to the the 
following minimization problem: 

Find the functions u,v which minimize the functional I(u,v) over all 
functions belonging to some admissible class where 

I(u,v) = // I[(pu) + (pv) - s] 2 + a [u - V + a)] 2 } dxdy 
ft x y u y x 

+ a j / (u - u Q ) 2 dy + a 2 / (v - v Q ) 2 dx, (9) 

and <Xq, and are some Lagrange multipliers. 
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Suppose that I(u,v) attains its minimum value for u = u* and 
v = v*. In Eq. (9), p, s and w are assumed to be known functions. If p 
is a positive function of x and y, then the original system (1) is 

elliptic. In the nonlinear case, when p depends on the solution u and 
v through Bernoulli's equation, the system is of mixed elliptic-hyperbolic 
type. 

We choose as our admissible class of functions those which satisfy the 

ft 

given boundary conditions and are twice dif ferentiable. Let u = u + 

^ a 

and v = v + tu, then 
■j[l(u,v) - I(u*,v*)] 


" l ! h Ve[<e u \, + - S J + “oKy - \y 

- JJ e 2 *2 iptcpu*)^ + (pv*) yy - s y ] + a 0 [v^ - u yx 

+ / e i "H ^ { P [ ( P u ) x + (pv ) - s] + a^(u - u Q )}dy 

+ I e 2 ^{pUp* 1 *^ + (pv *^ y ~ s ] + a 2 (v * " v 0 )l dx 

+ «0 / e i ^[“y - v x + “]dx 

- a o / e 2 " V x + ^ dy + 0(E 1 ) + 0(e 2 ) - 


+ w y ]}dxdy 
- w x ]}dxdy 


( 10 ) 


Since, by definition, vanishes along x = 0 and x = JL and r) 2 

vanishes along y = 0 and y = h, the third and fourth integrals in (10) are 
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identically zero. Following standard minimization arguments we can show that 
the functions u , v which minimize I(u,v) are the solutions of the 
boundary value problems shown in Figure 3 and 4, taking <Xq = 1. 


u = v - 0 ) 


y x 



u = v - 0) 

y x 


Figure 3. Pictorial Representation of the u Problem (Problem I) 



Figure 4. Pictorial Representation of the v Problem (Problem I) 




-7- 


Phillips [3] shows that when p — 1 the least squares formulation has a 
unique solution. Moreover, this solution satisfies the original problem. 

In an iterative procedure it is natural to use the problems defined in 
Figures 3 and 4 to solve for u and v, respectively. There are many ways in 
which iterative methods can be implemented to solve such a coupled system of 
equations. 

For Problem II, where the tangential rather than the normal component of 
the velocity is specified around the boundary, the functional we need to 
minimize takes the form 

J < u « v ) = // {[(pu) x + (pv) y - s] 2 + a 0 [u y - v x + w] 2 }dxdy 

+ “]_/(▼“ v Q ) 2 dy + a 2 / (u - u Q ) 2 dx. (11) 

In a similar fashion we can proceed to show that the functions u*, v* that 
minimize J(u,v) over all admissible functions must solve the boundary value 
problems, shown in Figures 5 and 6, again taking o Q = 1. 



(pu) x = -(pv) + s 


u = u r 


Figure 5. Pictorial Representation of the u Problem (Problem II) 



(P V )y = “(P U ) X + S 



(pv) y = “(P u ) x + S 


Figure 6. Pictorial Representation of the v Problem (Problem II) 


3. MODIFIED LEAST SQUARES METHOD 

If the original system is written in the form 


DU = B 


where 
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th e n the above classical least squares formulation leads to a system which can 
be represented symbolically by 

EDU = EB , (13) 


where 




To avoid the need for differentiation, a new variable V is introduced 
through the relationship 


Hence, Eq. (12) becomes 


DEV = B. 


(15) 


This choice is equivalent to applying a weighted least squares to the system 
(12) where the weight is related to 


R = (DE) -1 . 


The system of equations (15) is similar to the system (13) with one exception; 
the nonhomogeneous terms are not differentiated. With reference to Hafez [2], 

a positive weighting function Q may be introduced and Eq. (15) is slightly 
modified to read 
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DQEV = B. 


(16) 


In the particular example under consideration, Q is chosen to be 


Q = 




> 


and so V can be identified in terms of potential and stream functions. More 
precisely, we have 


v v + v p v 


DQE 


V 1/p V + a y (1/p V' 


Also, if V T = (c|),cp) , then Eq. (14) can be written in the following component 
form to express u and v in terms of <t> and 4> ; 


<!>x + 1/P = u 

4> y - 1/P «l» x = v. 


(17) 


For the case in which u> vanishes identically, the boundary conditions on 
<p and <l> can be chosen so that 4> vanishes identically and Eq. (16) reduces 

to 

(P ♦,), + (P *y)y 


= s. 


(18) 


- 11 - 


With the associated Neumann boundary conditions for Problem I or Dirichlet 
boundary conditions for Problem II, the formulation is complete. 

Alternatively, if s vanishes identically, the boundary conditions on <)> 
and i{> can be chosen so that <|> vanishes identically in £2 and Eq. (16) 
reduces to 

0f> /p) v + 0l> /p) v = -u> (19) 

x a y y 

and together with Dirichlet boundary conditions for Problem I or Neumann 
boundary conditions for Problem II, the formulation of this problem is 
complete. 

In the general case when neither s nor a) vanish identically in ft the 
governing equations are (18) and (19). We note that the problems for <f) and 
decouple except for the boundary conditions. 

This modified formulation is related to the Helmholtz theorem (see [2]) 
which allows a vector to be decomposed into two vectors, the first of which is 
curl free while the second is divergence free, i.e. , 

q_ = Vcf> + V x A. (20) 

For two-dimensional flows, the vector _A can be represented by one component, 
while for three-dimensional flows, at least two components of JV are needed. 

With the present method, s and w (or both) can be discontinuous. In 
such a case, unlike u and v, $ and ^ remain continuous. Their 
derivatives, of course, are not continuous. 

The jump conditions of a weak solution of Eqs. (18) and (19) are 
(assuming s and o> are integrable). 
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K1 - (£)„ [p* t ] - o (2D 

and 

l> x /p ] " ^ y / p ^ = °» (22) 

and since it is assumed that 4> = § and cj> = c|> , across a 

xy yx xy yx 

discontinuity D, the following relations hold 

[♦y] + t* x ] ‘ ° < 23 ) 

and 

[♦,] + (§)» [♦,] - °- w 

It is obvious then, that a linear combination of Eqs. (21), (22), (23), and 

(24) satisfy the jump conditions admitted by Eqs. (1), namely, 

[pu] - (§} D [pv] - 0 (25) 

and 

[v] + #) D [u] = °- < 26 ) 

In a sense, 4> and <l> are integrals of u and v and thus a second-order 
system can be constructed without differentiation. 


4. NUMERICAL RESULTS 

We present results for the case in which p - 1, i.e., Eqs. (1) and (2) 
are the Cauchy-Riemann equations. Equations (18) and (19) can be discretized 
using finite difference or finite element methods and the resulting algebraic 
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system may be solved using SOR, ADI, conjugate gradient or raultigrid 
techniques, for example* 

The system of Eqs* (18) and (19) represents two Poisson equations, one for 
$ an< 3 one for when p = 1. Many of the standard iterative methods may 
be used to yield a solution to these equations* We have chosen to use the 
multigrid method since such methods are fast and efficient for these problems. 
The domain Q is covered with a square grid of mesh length h = 1/N 
where N is a positive integer* Each of the Poisson equations: 


V 2 * = S, 

V 2 = -to , 


(27) 


is discretized using standard second-order central difference approximations. 
We restrict ourselves to Problem I in which the normal component of the 
velocity is specified around the boundary. We consider three types of problem 
and associated with each one will be different boundary conditions: 

(a) s = 0. In this case the boundary conditions can be chosen so that <f> 

vanishes identically. The resulting problem for is one in which 

Dirichlet conditions are given on the boundary. These are obtained by 
integrating the given velocity boundary conditions around . 

(b) w = 0. In this case the boundary conditions can be chosen so that 

vanishes identically. The resulting problem for <j> is one in which 

Neumann conditions are given on the boundary. Here we require that the 
compatibility condition be satisfied in order for a unique solution to 


exist. 
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(c) s £ 0, a) i 0. This is the general case where we need to solve for both 
(j) and i|j. The treatment of the boundary conditions in this case must 
be consistent with Eqs. (5) and (8). 


We consider a multigrid method of solution to these problems using the 
correction storage algorithm of Brandt [1]. Let Gj,»**,G m be a sequence of 
grids approximating the domain ft with corresponding mesh sizes 
Let h k = 2h k+1 for k = l,»»»,m-l. The problem is discretized on each grid 
G k as described above. We use the same components in the multigrid procedure 

as those chosen by Phillips [3]. Very briefly these are pointwise Gauss- 

Seidel with the points ordered in the checkerboard manner for relaxation, 

half-weighting to transfer the residuals to the coarser grid and bilinear 
interpolation to transfer the correction to the fine grid. 

We consider three test problems defined in the unit square and 

characterized as follows: 


(i) s = 0, 


o) 


TT 2 2 

= f— ] (m + n ) sin(rtJTx) sin(mny), 
v m' 


u = sin(rorx) cos(mny), v = cos(tmx) sin(miTy), 


♦ = o, 


to = -f— 1 sin(mix) sin(mny); 
^ nnr ' 


(ii) 



cos(nnx) cos(mrry). 


a) = 0, 


u = sin(mrx) cos(miry), 


v = f— ) cos(nirx) sin(nnry), 
v Tr 
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<t> = cos (mix) cos(nncy). 


<P = 0; 


(Hi) s = (m+n)x cos (nnx) cos(mxy), to = (m-n)it sin(nxx) sin(nmy). 


u - sin(nxx) cos(mity), v = cos (nnx) sin(nmy), 


4> = — 2 + 2 cos (nnx) cos(mny), <p = ■ sin(nnx) sin(mny), m * n. 
(m +n )n (m +n )n 

In all of the numerical experiments we took the finest grid to be such 
that N = 32 and considered a total of five grids in the multigrid context. 
The initial approximation was taken to be zero everywhere except where the 
solution was specified by the boundary conditions. The iterations were 
terminated when the i^-norm of the residual had been reduced by a factor of 
10 from its initial value. After solving Eq. (27 ) for 4> and/or cp, we 
obtain the corresponding values of u and v by approximating the left-hand- 
side of Eq. (20) by central differences. We note that increased accuracy 
could have been obtained if we had used some global interpolation scheme. 

The first problem corresponds to one of the examples considered by 
Phillips [3]. Results very similar to those for u and v discussed in [3] 
are obtained using the modified least squares approach outlined in this paper. 
Close agreement was achieved in measurements of different norms of the errors 
in the discrete solution and in the convergence histories of the two 
techniques. Similar conclusions were reached for the second problem. 


-16- 


With the normal component of the velocity specified around the boundary, 
the boundary conditions in terms of <t> and <J> take the form: 


4> x + 4> y = u Q along x = 0, x = 1, 

<t> y “ 4» x s v Q along y = 0, y = 1. 


(28) 


We have chosen to decouple these boundary conditions in terms of <1> and 
by putting 4> = 0 on the boundary and obtaining conditions for <J> from Eq. 
(28). 

The details of the algorithm applied to the third problem are given in 
Tables I, II, and III. In Table I we give the details of the <|> calculation. 
We give the number of work units required to attain the convergence criterion, 
and the asymptotic convergence factor X for various values of m and n. 

Similar information is furnished in Table II for the i|> calculation. In 

Table XII we give norms of the errors in the velocities u and v for 
various values of m and n. We have used the notation a.b — c for 
a.b x 10 c . 

It can be seen that the method exhibits the usual multigrid behaviour by 
examining the asymptotic convergence factors obtained. The accuracy of the 
discrete approximation to the third problem decreases as m and n increase 
as one might expect since the number of mesh points per wavelength of the 


solution decreases. 
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Table I. Details of the <j> Calculation 


m 

n 

wu 

X 

i 

2 

21.41 

0.56 

2 

1 

21.41 

0.56 

3 

1 

17.88 

0.52 

1 

5 

17.25 

0.49 


Table IX. Details of the (J; Calculation 


in 

n 

WU 

X 

1 

2 

18.13 

0.55 

2 

1 

18.13 

0.55 

3 

1 

17.81 

0.52 

1 

5 

17.50 

0.49 


Table III. Error Norms of the Solution 


nun 

oo 

lulj. 

Hull 2 

11 v 11 

CO 

HvH x 

" vll 2 

0.46-2 

0.18-2 

0.23-2 

0.81-3 

0.31-3 

0.39-3 

0.46-2 

0.18-2 

0.22-2 

0.81-3 

0.31-3 

0.39-3 

0.10-1 

0.40-2 

0.50-2 

0.27-2 

0.11-2 

0.13-2 

0.12-1 

0.46-2 

0.58-2 

0.27-1 

0.10-1 

0.13-1 


1 
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5. CONCLUDING REMARKS 

An application of the least squares method to a system of first-order 
equations has led to a formulation in terms of second-order equations. The 
process has also determined the additional boundary conditions necessary for 
the higher order equations. The duality of problems in which either the 
normal or the tangential components of velocity are specified on the boundary 
has been indicated. Also, a modified least squares approach has been 
outlined. This approach avoids the need to differentiate the source or 
vorticity functions. Numerical examples demonstrate the effectiveness of the 
method. 

Future work will concentrate on treatment of the nonlinear problem, i.e., 
when p = p(u). Possible extensions of this technique to three dimensions 


using either two or three stream functions are under examination. 
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